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We provide a recursive method for constructing product formula approximations to exponentials 
of commutators, giving the first approximations that are accurate to arbitrarily high order. Using 
these formulas, we show how to approximate unitary exponentials of (possibly nested) commutators 
using exponentials of the elementary operators, and we upper bound the number of elementary ex- 
ponentials needed to implement the desired operation within a given error tolerance. By presenting 
an algorithm for quantum search using evolution according to a commutator, we show that the scal- 
ing of the number of exponentials in our product formulas with the evolution time is nearly optimal. 
Finally, we discuss applications of our product formulas to quantum control and to implementing 
anticommutators, providing new methods for simulating many-body interaction Hamiltonians. 



I. INTRODUCTION 

Product formulas provide a way of approximating a single operator exponential with a product of simpler operator 
exponentials. Such formulas are useful in numerical analysis, where they can be applied to the solution of differential 
equations (see for example [1]). More recently, product formulas have become a key tool for quantum simulation 
[2-10]. Hamiltonian simulation using product formulas has numerous applications in quantum information process- 
ing, including simulating quantum mechanics [ ], implementing continuous-time quantum algorithms [4, 11-13], and 
quantum control techniques (see for example [14]). 

The primary application of product formulas is to represent exponentials of sums. Although exponentials of commu- 
tators are not as ubiquitous, they arise naturally via their role in Lie groups. Exponentials of commutators appear in 
numerous asymptotic expansions, including the Baker-Campbell-Hausdorff series and the Magnus expansion. They 
also play a role in quantum computation, such as in the Solovay-Kitaev theorem [15], which constructively proves 
that any finite universal gate set is sufficient to perform efficient universal quantum computation, and in quantum 
control [14], where product formulas for exponentials of commutators can be used to suppress couplings or introduce 
ones that are not naturally present. 

Although the theory of product formula approximations for exponentials of sums is well understood, it is con- 
siderably less developed in the case of commutator exponentials. Product formula approximations to commutator 
exponentials approximate an exponential of the form exp([A, B]T) for operators A and B and a real number T with a 
sequence containing exponentials of A and of B. In the limit of small T, low-order product formulas for exponentials 
of commutators are well known, but we are not aware of any previous systematic method of constructing high-order 
approximations. 

In this work, we construct arbitrarily high-order product formula approximations to exponentials of commutators. 
Our formulas are analogous to Suzuki's seminal work on product formulas [16], but apply to the case where the 
exponentiated operator is a commutator, rather than a sum, of two operators. We do not explicitly consider cases 
where the operator is a sum of commutators or is an ordered commutator exponential, but such cases can be addressed 
by combining our results with existing product formula approximations for exponentials of sums [16] (see [8] for an 
improved analysis in some cases) or ordered operator exponentials [7, 9, 17]. 

Our results provide a method to simulate exponentials of the form e^ ,s l T , for any desired T € K, using devices 
that can enact evolution under A or B separately. Specifically, we imagine that we have a pair of devices T>a and T>g 
that take as input an evolution time t and perform operations e At and e Bt , respectively. Physically, we can imagine 
that T>a and T>b represent control fields that enact a desired evolution in a quantum system and t represents the 
time for which those control fields are applied. Alternatively, we can imagine that iA and iB are Hamiltonians that 
can be easily simulated and that T>a and T>g represent quantum simulation algorithms performing the corresponding 
evolution. In either case, we measure the efficiency of our formulas by the number of times these devices need to be 
used to simulate e^ A ' B ^ T . Note that our figure of merit is not the total amount of time the control fields are applied: 
an elementary evolution e At or e Bt has unit cost independent of t. 

The remainder of the paper is organized as follows. In Section II we present basic product formulas for approximating 
exponentials of commutators. We combine these formulas to approximate exponentials of nested commutators in 
Section III. Upper bounds for the approximation errors are derived in Section IV. These error bounds are applied in 
Section V, where we show that the resulting formulas use a number of exponentials that is only slightly superlinear 
in the evolution time. We show that this performance is nearly optimal in Section VI by proving that sublinear 
simulation would violate the quantum lower bound on the query complexity of unstructured search. We present 
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applications of our techniques in Section VII, including simple examples of quantum control as well as a method for 
simulating exponentials of anticommutators of operators. In particular, the simulation of anticommutators provides 
a novel method for implementing many-body interactions in quantum systems. We conclude in Section VIII with a 
discussion of the results and some open problems. 



II. BASIC PRODUCT FORMULAS 

In this section we present basic formulas approximating the operator e^ A ' B ^ + for small t as a product of powers 
of e At and e Bt , for a given positive integer k. Choosing k = 1 yields the most efficient formulas for the case where 
B can be implemented directly. We consider higher values of k for the case where B is itself a (possibly nested) 
commutator, as discussed further in Section III. We present two recursive constructions, one for k odd and another 
for k even, giving high-order approximation formulas for e^ A ' B ^ in terms of exponentials of A and B. For every 
integer p > 1, we present a formula with approximation error 0(t 2p+k+1 ) in the limit of small t. 

Since the k — 1 case is the most natural, we begin with splitting formulas for the case where k is odd. We then 
discuss the simpler case where k is even in Section II B. 



A. Odd-fc Formulas 



We now develop a recursive approximation-building method that can be used to construct an arbitrarily high-order 
approximation to e^ A ' B ^ tk+1 in terms of a product of powers of e At and e Bt , where k is odd. The construction uses 
an initial approximation to the time evolution, as follows. 

Lemma 1. Let A and B be bounded operators, let k > 1 be a real number, and define 

V 1<k (At, Bt k ) := e M e Bt " e- At e- Bt " . (1) 

Then V hk (At,Bt k ) = e {A,B]t^+0{t*+>) _ 

Proof. The Baker-Campbell-Hausdorff (BCH) formula implies 

Vi k (At, Bt k ) = e At+Bt k + ±[A,B]t k + 1 +0(t k + 2 ) e -At-Bt k + ±iA,B]t k + 1 +0(t k + 2 )^ 

A second use of the BCH formula gives 

V hk (At,Bt k ) = e [^.s]t fc+1 +o(t fe+2 ) ! (3) 
which completes the proof. □ 

The product formula V]_ t i(At,Bt) = e At e Bt e~ At e~ Bt is known as the group commutator. This formula has many 
applications, including generating optimal control sequences [14] and approximating unitary gates via the Solovay- 
Kitaev theorem [15]. We show that higher-order generalizations of this formula can be constructed using an iterative 
approximation-building method that is reminiscent of Suzuki's method [16]. 

To use our technique, we must have product formulas for the inverses of our approximations. Our approximations 
built from V± lk (At, Bt k ) possess one of two symmetry properties that make their inverses simple to compute. 

Definition 1. A product formula U is symmetric if U(X,Y) = U(Y,X)^ 1 and is antisymmetric if U(X,Y) = 
U(—Y, — X)^ 1 , for all bounded operators X and Y . 

In particular, V\ ik is symmetric, and we will see that high-order approximations constructed using this product 
formula are either symmetric or antisymmetric. 

Now we are ready to describe the main result of this section, which shows how to construct an arbitrarily high-order 
approximation. 

Theorem 2. Let A and B be bounded operators, let k > 1 be an odd integer, and let V Ptk (At, Bt k ) be a product 
formula with V p>k (At, Bt k ) = e^ A,B ^ + +°(* P+ ) for some positive integer p. Let 

V p+hk (At,Bt k ) := V p , k {A lp t,B{ lp t) k )V p , k {-A lp t,-B{ lp t) k ) 

x V Ptk (Ap p t, Bi^tfy^i-Appt, —B(/3 p t) k )~ 1 

x V pM {A lp t, B{ lp t) k )V p , k {-A lp t, -B{ lp t) k ) (4) 
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where 



^:=(2r p )V(W), 7p:=(1/4 + rp) V(W), r p :^ v 

4(2- 2 2 p+ fc + 1 



T/ien Vp +M (Ai,.Bi) = exp(L4, B]t 2 + 0{t 2 ^+^ +1 )) and V p+1 , k (At, Bt k ) = exp([A, B]t k+1 + 0(< 2 (p+ 1 )+ fe + 1 )) j/jfe > 1. 
Furthermore, V p ^i -k is symmetric ifV Ptk is antisymmetric and is antisymmetric ifV p ^k is symmetric. 

Proof. The assumption that \\U(t) - V p<k (At, Bt k )\\ € 0(t 2p+k ), where U(t) := e [A t B]t k+1 ^ i mp ii es that there exist 
operators C{A,B), D(A,B), and E(A,B) such that 

V p ,k(M, Bt k ) = U(t) + C(A, B)t 2p+k + D(A, B)t 2p+k+1 + E(A, B)t 2p+k+2 + 0{t 2p+k+3 ). (6) 

Every term in the Taylor series of V p ^(At, Bt k ) is a product of powers of At and Bt k , each of which contributes 
an odd power of t. Thus each term in C(A,B) and E(A,B) contains an odd total number of A and B operators. 
Similarly, because 2p + k + 1 is even, each term in D(A, B) contains an even number of A and B operators. Therefore, 

C(A,B) = -C(-A,-B) 
D(A, B) = D(-A,-B) 

E(A,B) = -E(-A,-B). (7) 

We use these properties to simplify V p>k (At, Bt k )V P: k{~ At, — Bt k ) for arbitrary t. Specifically, we expand V p ,k as a 
power series and use (6) to show that 

V Ptk (At,Bt k )V p , k (-At,-Bt k ) = U(t) 2 + [C(A,B)U(t) - U{t)C{A, B)]t 2p+k 

+ [D(A, B)U(t) + U(t)D(A, B)]t 2p+k+1 

+ [E(A, B)U{t) - U{t)E(A, B)]t 2p+k+2 + 0{t 2p+k+3 ). (8) 

We then Taylor expand each U(t) (but not U(t) 2 = U(2 1 ^ 1+k h)) in (8) to lowest order in t and find that 

V Ptk (M, Bt k )V pM {-At, -Bt k ) = U(t) 2 + 2D(A, B)t 2p+k+1 + 0{t 2p+k+3 ) + 0{t 2p+2k+1 ). (9) 

This implies that (4) has no error terms of order t 2p+k because it is a product of three pairs of product formula 
approximations of the same form as (9). 

Next we show that the careful choice of [3 P and 7 P eliminates the term of order t 2p+k+1 . To do so, we relate the 
error terms of V p ^ k (At, Bt k ) and V p ^(At, Bt k )~ 1 . Similarly to (6), we have 

V Plk (At, Bt k y l = U{t)- 1 + C(A, B)t 2p+k + D(A, B)t 2p+k+l + E(A, B)t 2p+k+2 + 0(t 2p+k+3 ) (10) 

for some operators C(A,B), D(A,B), and E{A,B). This expansion directly follows from the symmetry properties 
of Vp fc. For example, if V p ,% is symmetric, then V pk {At, Bt k )~ 1 = V p ^ k (Bt k , At). We know from the previous 
discussion that e [B ' A]tk+1 = V Pik {Bt k ,At) + 0{t 2p+k+1 ) = E/(i) _1 . The symmetry of the formula then implies that 
C/(<) _1 = V p , k (At, Bt k y x + 0(t 2p+k+1 ), which justifies (10). The anti-symmetric case follows similarly. 

Since each term in C(A, B) and E(A, B) consists of an odd number of A and B operators and each term in D(A, B) 
consists of an even number of such operators, we have (similarly to (7)) 

C(A,B) = -C(-A,-B) 
D(A, B) = D(-A,-B) 

E(A,B) = -E(-A,-B). (11) 

Equations (4), (9), and (10) then imply that 

V p+1 MAt,Bt k ) = (U{ lp t) 2 + 2D{A,B){ lp t) 2p+k+1 ) 

x (U{fi p ty 2 + 2D(A,B)(/3 p t) 2p+k+1 ) 

X (U{ lp t) 2 +2D{A,B)( lp tf P+k+1 ) + 0(^in{2p+fc+3,2p+2fc+l}) 

= U(t) + U~f 2p+k+1 D{A, B) + 2(3 2p+k+1 D(A, B)] t 2p+k+1 + o(t mini - 2p+k+3 ' 2p+2k+1 ^). (12) 
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We can relate D(A, B) to D(A, B) by noting that 

1 = V p , k (At,Bt k )V Ptk (At,Bt k )- 1 
= (U(t) + C(A, B)t 2p+k + D(A, B)t 2p+k+1 + E(A, B)t 2p+k+2 ) 

x (Uity 1 + C(A, B)t 2p+k + D(A, B)t 2p+k+1 + E(A, B)t 2p+k+2 ) + 0(t min ^ p+k+3 ' 2p+2k+1 ^). (13) 

By Taylor expanding the resulting formula, we find D(A,B) = —D(A,B) (as well as C(A, B) = —C(A,B), and if 
k > 1 then E(A, B) = —E(A,B)). 

Recalling the definitions of f3 p and 7 p , we can substitute D(A, B) = —D(A, B) into (12) to find 

V p+hk {At,Bt k ) = U(t) + Uil/A + rp) 2 - 1 ^ 1 -2(2r p ) 2 - I ^r 1 } D{A, B)t 2p+k+1 + 0(t min V p+k+3 > 2p+2k+1 ~>). (14) 

The value of r p specified in (5) is a root of the expression in square brackets, so the term of order £ 2p+fc+1 vanishes. 
This demonstrates that 

V p+hk (At, Bt k ) = e [^]' fc + 1 + ( t nxin{2 P +fc+3,2 P +2fc+l} ) (15) 

as claimed. 

Finally, we must show that V p+ \_ k is either symmetric or antisymmetric. First, suppose V Ptk is symmetric. Then 
V p+1 . k {At, Bt k )V p+hk (-Bt k , -At) 

= (^fe(^7 P i,^(7 P fe )^,fc(-v4 7p i,^B( 7p t) fc )^ fc (B(/? p i) fc ^/3pi) 
x V pM {-B{p p t) k ,-Ap p t)V p , k {A lp t,B{ lp t) k )V p , k {^ 
x (V p , k (-B( lp t) k , -A lp t)V p , k (B( lp t) k ,A lp t)V p , k (-A/3 p t, -B(j3 p t) k ) 

x ^, fc (A/3 p ^,s(/3 p ^) fe )l/ p , fc (-B(7 P ^)^-A7 P ^)^ fc (B(7 P ^)^A7 P ^)) = i. (16) 

Thus Vp+i^ is antisymmetric. A similar calculation shows that if V Ptk is antisymmetric then V p +\ yk is symmetric. □ 

The result of Theorem 2 shows how to recursively construct the product formula V p . k starting from Vt ik - Note that 
this construction involves terms of the form V q ^ k (AXt, B(\t) k )~ 1 . Since V q , k is either symmetric or antisymmetric, its 
inverse can be represented explicitly using Definition 1. 

The following corollary shows that symmetrization alone can increase the order of the approximation of V p _\ from 
0(t 2p+1 ) to 0(t 2p+2 ) at the cost of doubling the number of exponentials. 

Corollary 3. Let A and B be bounded operators, let V p ,i(At, Bt) satisfy V p ,i(At, Bt) = e^ A ' B ^ P+ and define 

V p l (At, Bt) := V pA (At/V2,Bt/V2)V p ,i(-At/V2, -Bt/V2). (17) 

Then 

V; s (At,Bt) = e lA ' B]t2+ ° (t2P+2) - (18) 

Proof. This is a simple consequence of (9) with k = 1. Since every term of order 2p+ 1 in V p _i (At, Bt) must contain an 
odd total number of As and Bs, the terms of order 2p+ 1 in V Py i(At, Bt) and V Pt i(— At, —Bt) are equal and opposite. 
By Taylor expanding V p ^(Atj\j2, Btj \f2)V p ^(— Atj\[2, —Bt/y/2), it is easy to see that the error terms of order 2p+ 1 
cancel. The error in V pl (At,Bt) is therefore 0(t 2p+2 ) as claimed. □ 

Note that a similar symmetrization would improve the error bound of V P , k from 0(t 2p+k ) to 0(t 2p+k+1 ), but this 
gives no improvement over Theorem 2 for k > 1 since that theorem already shows that the error is 0(t 2p+k+1 ). 

For k > 1, we can apply Theorem 2 recursively to produce a formula V p . k with 4 x 6 P_1 exponentials having error 
0(t 2p+k+1 ). For k — 1, applying Theorem 2 recursively followed by one application of Corollary 3 gives a formula 
V pl with 8 x 6 P_1 exponentials having error 0(t 2p+2 ). 

Figure 1 presents a numerical example showing improved error scaling as p increases. This example considers 
simulating the commutator of the Pauli operators a x and a z using k = 1. These data suggest that our upper bound 
on the error scaling is in fact tight. 
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FIG. 1. Error scaling of the formulas V Pl i(— ia x t, —ia z t) for p = 1, 2, 3 as approximations of e ~^ m '' 7 ^ t , where || • || is the 2-norm. 



B. Even-fc Formulas 

If k is even then the above approach does not apply, so different reasoning must be used to generate high-order 
approximations to e^ A ' B ^ f . We show, somewhat surprisingly, that Suzuki's recursive approximation-building for- 
mula [ ] for the exponential of a sum can be used to generate arbitrarily high-order approximation formulas for 
commutators from a relatively simple initial approximation. The initial formula is as follows. 

Lemma 4. Let A and B be bounded operators, let £ k := 2^ 1 ^ 1+k \ let k > be even, and define 

W 1:k (At,Bt k ) := e At ^e Btk ^e- 2At ^e- Btk ^e At ^. (19) 

Then 

e [A,B]t k+1 = w hk (At, Bt k ) + 0(t k+3 ). (20) 

Proof. Since U(t) := e^- A ' B ^ tk+1 = e^- B '~ A ^ tk+1 , it follows from Lemma 1 (which applies for any k) that 

g[A,B]t k+1 = VlJc{A t^ kl Bt k e k )V hk (Bt k e k ,-At^) + 0(t k+2 ) 

= (e^eBt^V^e- 3 '^) (V^V^e-^e^*) + 0(t k + 2 ). (21) 

The dominant term of equation (21) reduces to Wi, k (At, Bt k ) by simplifying the middle terms. 

It remains to show that the term proportional to t k+2 vanishes. This follows from the fact that W\ t k[At, Bt k ) = 
W l<k (-At, Bt k )- X . Thus, defining C(A, B) via W hk (At, Bt k ) = U(t) + C(A, B)t k + 2 + 0{t k+3 ), we have 

1 = W hk (At, Bt k )W hk {-At, Bt k ) = 1 + (U{t)C{-A, B) + C(A, B)U{-t))t k+2 + 0(t k+3 ). (22) 

Because k is even, C must be composed of an even number of As, so C{— A,B) = C(A,B). Since (22) must hold for 
arbitrary t, we have C(A, B) — 0. □ 

The product formula in Lemma 4 is reminiscent of the Strang splitting formula, and indeed reduces to that approx- 
imation to 1 = e^ A_j4 ^ when B = 1, or if we take B — B k for some operator B and substitute k = into the resulting 
formulas. Since Suzuki's iterative approximation-building method can refine the Strang splitting into arbitrarily high- 
order formulas, one might suppose that the same recursion could approximate exponentials of commutators. We make 
this intuition precise in the following theorem. 
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Theorem 5. Let A and B be bounded operators, let k > be an even integer, let p > 1 be an integer. Suppose that 
W Ptk {At,Bt k ) = W p ^ k (-At,Bt k )- 1 and W p , k (At,Bt k ) = e [AB]* fc+1 +0(t 2p+fc+1 ) f and define 

W p+hk (At, Bt k ) := W p>fc (AV,S(V) fc ) 2 %,it(-^M P ^-B( Mp t) fc )W p , fe (AV,S(i/ p t) fc ) 2 (23) 

where 

fc+i 

M P := (4 Sp ) 1/(fc+1) , v p := (1/4 + s p )^ k+i \ s p := (24) 

4(4 — 4 2 p+ fc + 1 ) 

T/ien 

W p+hk (At, Bt k ) = exp ([A, B]t k+1 + 0(t 2 P+ k + 3 )) . (25) 

Proof. We follow the same reasoning used to analyze the Suzuki formulas. We have 

W p ,k{At, Bt k ) = U(t) + C(A, B)t k+2p+1 + 0(t 2p+k+2 ) (26) 

for some C(A, B). Substituting (26) into (23) gives 

W v+l<k (At, Bt k ) = U(t)+ (4C(A, B)is 2p+k+1 + C(-A, B)nf +k+1 ) t 2p+k+1 + 0{t 2p+k+2 ). (27) 

Each term comprising C(A, B) must contain an odd number of As, because each A is associated with t and each B 
is associated with t k , so since k is even, a term proportional to t 2p+k+1 can only be formed from an odd number of 
As. Thus C(A,B) = -C(-A,B). Therefore, the coefficient of C(A, B) is zero when 

4(1/4 + Sp )(2p+fc+i)/(fc+D _ (4Sp) (2 P +fe+i)/(fc+i) = a (28) 

The value of s p specified in (24) is a root of (28) and hence 

\\U(t) - W p+l , k {At,Bt k )\\ G 0(t 2p+k+2 ). (29) 

Finally, we show that by symmetry properties of W p +\ :k , the actual error scaling is better than in (29). We have 

W p+ i ik (At, Bt k ) = U(t) + D(A, B)t 2p+k+2 + 0(t 2p+k+3 ). (30) 

Since W Pik {At, Bt k ) = W p . k (-At, Bt k y x , it is easy to see by multiplication and (23) that W p+ i ik {At, Bt k ) = 
W p+ i t k(-At,Bt k )~ 1 . We conclude that D(A, B) = by the same argument used in (22). Since D(A, B) = 0, 
we have W p+hk {At, Bt k ) = U{t) + 0(t 2 ^ p+1 '>+ k+1 ) as claimed. □ 

The output from the recursive formula in Theorem 5 can be used as input, so the theorem shows how to refine W± k 
into an arbitrarily high-order approximation. The basic formula Wx tk has 5 exponentials, and each iteration increases 
the number of exponentials by a factor of 5, so W Pik {At, Bt k ) consists of 5 P exponentials. 

Figure 2 presents a numerical example showing improved error scaling as p increases with k = 2. As in the case of 
odd k, these data suggest that our upper bound on the error is tight. 

III. EXPONENTIALS OF NESTED COMMUTATORS 

In general, we may be interested in approximating exponentials of nested commutators of the form 

Z k := [Ak, [Ak-i, [■ • ■ , [Ai, Aq] . . .]]] (31) 

(i.e., Zq — Aq and Z k = [A k ,Z k _i] for k > 0). Our strategy for approximating such exponentials is simple. For 
example, suppose k is odd. We first approximate e Zkt with V P}k (A k t, Z k ^\t k ). The resulting expression still 
contains commutator exponentials because of the presence of exponentials of Z k _\. We approximate each exponential 
of the form e Zk - 1 ^ t ^ by W p _ k ^i{A k _i\t, Z k _ 2 (^t) k ~ 1 ), leaving exponentials of Z k ^ 2 in the formula. This process 
is repeated recursively until no exponentials of commutators remain. The case where k is even can be addressed 
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similarly. In both cases, the resulting product formula is denoted lA p {Akt, . . . , Aot). In other words, we approximate 
a single commutator with U p (Ait, Aot) := V p i(Ait, Aot), and we recursively define U p for q > 1 as 

- I Z} • <*> 

for all 9 > 1, where Rep(x, a — > b) indicates that we replace every instance of a in x with 6. For example, 

Ui(A 2 t,A 1 t,A t) = Rep(W 1 , 2 (A 2 t,Z 1 t 2 ), e Zl(At)2 ->■ Wi(AiAi, AqM)) 

= Rep(e A2 « lt e Zl (« lt ) 2 e 2A2 « lt e- Zl ^ t ) 2 e A2 « lt , e Zl ( At ) 2 -> Wi(^iAi, A A<)) 

= e^Hi^fct, A eii)e 2A2?1 U(^i6t, AoCii)"^^ 1 *. (33) 

At first glance, it may seem surprising that a recursive approximation to e Zkt + can use approximations with error 
of order less than 2p + k + 1 (rather than using approximations with error 0(t 2p+k+1 ) at each stage in the recursion) 
when the desired overall error is 0{t 2p+k+1 ). However, since a term with relatively large error is multiplied by a 
relatively small term, this recursive process straightforwardly yields formulas with error 0(t 2p+k+1 ). 

Lemma 6. For any positive integer k, let Aq, . . . ,Ak be bounded operators, define Zk by equation (31), and let 
Up(Akt, . . . , A^t) be the product formula defined above. Then U p {Akt, . . . , Aot) = e Zht + +°( t F+ + ) . 

Proof. We use induction on k. By Corollary 3, U p (Ait,A t) — e z i t2 +°( t2p+2 ) ; establishing the base case. For the 
induction step, Theorem 2 and Theorem 5 imply that 

U p (A q t, A Q t) = exp([A q t, Z g _xt q + 0{t 2p+ i)]) 

= cxp(Z q t q+1 + 0{t 2p+q+1 )). (34) 

The desired result then follows by induction. □ 

To calculate the number of exponentials N Pt k appearing in U p (Akt, . . . , Aot), recall from Section II that V p 1 uses 
N p ,i — 8 x 6 P_1 exponentials. If k > 1 is odd, then V Py k uses 4 x 6 P_1 exponentials, half of which require further 



expansion into N p ,k-i exponentials. If fc is even, then W P: k uses 5 P exponentials, 2/5 of which require further expansion 
into N Py k-i exponentials. Thus, for fc > 1, we have 

AT _ /5P- 1 (3 + 2iV M _ 1 ) keven 

| 2x6P -i (1 + 7Vpfe i) fcodd W 

While it is cumbersome to present the solution in closed form, it is easy to see that for fixed k we have N Pj k — 0(6 pk ) 
since we increase the number of exponentials by a factor of 0(6 P ) with each iteration. 

The approach described above uses both odd-fc and even-fc product formulas. It is possible to derive similar formulas 
using only odd-fc product formulas, but their error bounds are more difficult to prove and the resulting approximations 
are less efficient. 



IV. ERROR BOUNDS FOR NESTED COMMUTATORS 



So far, we have presented product formulas that approximate the evolution according to a (nested) commutator to 
arbitrarily high order. We now provide simple upper bounds on the approximation error that results when using these 
formulas to simulate evolution for a sufficiently small time. We use the following notation to denote the remainder of 
a Taylor series expansion truncated at order v — 1. 

Definition 2. Let f{t) be a bounded operator for any tgK with a Taylor series f(t) — X^^Lo a nt n ■ Then 

oo 

!*.„(/(*)):= X>«*" ( 36 ) 

n— v 

We use the following standard result (with a proof included for completeness). 

Lemma 7. Let {A q : q = 1, . . . , N} be a set of bounded operators and let \\-\\ be a submultiplicative norm. Then for 
any positive integer v , 



r „ n e 

\<Z=1 

Proof. Since each A„ is bounded, Taylor's theorem implies that 



N 



N 



{A q ty 



(37) 



JV 



\g=l p=0 

The submultiplicativity of ||-|| and the triangle inequality imply that 

/ N 

<Rjn 

\q=l 



(38) 



\q=lp=0 



P 



E 



p=0 



N oo 



< 



MIIE 



(U q \\ty 



\q=l p=0 



as claimed. 



(39) 
□ 



The following lemma establishes general error bounds for exponentials of (possibly nested) commutators. 
Lemma 8. Let {Aj : j = 0, . . . ,k} be a set of bounded operators and define by equation (31). Suppose that 

1. U p {A k t, . . . ,A t) = e z J tk+1 +o(t") j or some integer v > k+ 1, 

2. U p (A k t, . . .,A t) = l\qZ! k e x * A ^ where N~l |A,| < Q P , k , 

3. A > 2||A,-|| for all j = 0, 



, k, 
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4- At < 



ln(2) 



N p , k Q p , k > 

5- N Pi kQ Ptk > 1. 
Then 



and 



e Zkt -U p (A k t,...,Aot) 



< 



Proof. By assumption 1 and the triangle inequality, we have 



U p {A k t,...,A t)\\ <||R„(e Zfct )\\ + \\TL„{U p {A k t,...,A t))\\. 



The proof of the lemma follows from upper bounds for both of these terms. 

We first bound WR, u (U p (A k t, . . . ,^4qO)II- Assumptions 2 and 3, together with Lemma 7, imply that 



||R,(e 



\ 1 A 1 t\ 2 A 2 t 



)|| < R„ (eM* MAJ* ) < R^e^.^M*) 



(40) 



(41) 



(42) 



so Taylor's theorem gives 



R 



< 



2(N Ptk Q Ptk Atf 



where we have used assumption 4 to simplify the bound. A standard variant of Stirling's formula, namely [18] 



n! > v27rnn"e 



gives 



\R v (U p {A k t,...,Aot))\\ < 



2 f eN n . k Q n . k kt\ v 



'2ttv 



V 



< 



The corresponding bound for the remainder term of e Zkt follows similarly: 

|R, (e^ +1 ) 

Since ||Z fc ||i fc+1 < (At) k+1 < In 2 (by assumptions 3, 4, and 5), we have 



( || Zfc ||i/(* + D f) (* + i)W ||Zfc||tfc+1 

r— t 6 
i fc+i i • 



(43) 



(44) 



(45) 



(46) 



R 



„ ( e ^ tfc+1 ) 



< 



< 



i fc+i i- 

2 



e i/(*+i)Af 



^/(fc + i) VW(fc + i)) 1/(fc+1) 



(47) 



where the second step uses (44) and At < 1 (by assumptions 4 and 5). To put this in a similar form to (45), we use 
assumption 5 to find 



R, 



< 



y/2nv/(k + 1) 



fc + 1 W(/c+i) / eNpkQpkAt y 



Since (fc + l)/e < 1 for > 1 and v/{k + 1) > 1 by assumption 1, we have 



so 



R„ 



(48) 



(49) 



(50) 



10 



For k > 1 we have (k + lf/ 2 e - k < 2 3 / 2 /e < 1.05, 



4.1 { eN p , k Q Pik At y 



^ V ^ 1 /( fe+1 ) 

where the second step uses assumption 5. Since 4.1/\/27r^ < 4.1/v / 67r < 1, the result follows. □ 

The error bound of Lemma 8 can be used to estimate the maximum size of a time step for product formulas given 
by Theorem 2 for k = 1 or any of the product formulas considered in Lemma 6, given upper bounds for N Pjk Q P ,k- 
Recall from the previous section that for constant k, N Pt k — 0(6 pk ). Furthermore, since r p and s p are monotonically 
decreasing functions of p (for fixed fc), it can easily be seen that they scale as 1/2 + o(l) and 1/4 + o(l), respectively. 
This similarly implies that j3 p = lp = (l/2) 1 /( fe+1 ) + o(l) and fi p = v p = (l/3) 1 /( fe+1 ' + o(l) for fixed k. Since Q p , fe 
is the maximum value of products of such terms, it follows that Q p ^ k & 0(1) and thus N p ^Q P ,k & 0(6 pk ). In fact, a 
more detailed analysis shows that N Ptk Q p , k & 0((2v / 6) p ' c ), but this does not significantly improve our results. 

Note that the bound of Lemma 8 does not directly apply to product formulas that result from using k > 1 in 
Theorem 2, or any of the formulas considered in Theorem 5, because such product formulas contain varying powers 
of t. Error bounds for such cases could be proved by similar reasoning, but are unlikely to give bounds that are as 
tight, and the cases covered by Lemma 8 are more relevant for our applications. 

V. NEAR-LINEAR SCALING WITH TOTAL EVOLUTION TIME 

Lemma 8 bounds the distance between the true evolution U(t) and the approximation U p {A k t, . . . , A t). Although 
the formulas are very accurate for small t, the error (as measured by a sub-multiplicative norm such as the spectral 
norm) can be large if At > 1. To construct a product formula approximation that is close to U(t) for general At, we 
divide the evolution into short segments and show that the cumulative error from these segments is at most e. 

To ensure additivity of approximation errors from each time step, we restrict ourselves to the case where the 
operators Aj are anti-Hermitian (as in quantum mechanics), a restriction that was unnecessary in previous sections. 
By carefully choosing p to be a function of t and e, the number of exponentials needed to approximate U(t) within 
error e scales only slightly superlinearly with t k+1 (the total evolution time) and sub-polynomially with 1/e. 

The following lemma is the main result of this section. It provides an upper bound for the minimum number of 
time steps, r, needed to ensure that the cumulative approximation errors from all r time steps sum to at most e. 

Lemma 9. Let {Aj : j = 1, . . . , k} be a set of bounded anti-Hermitian operators and let e > 0. // Assumptions 3, 1, 
and 2 of Lemma 8 are satisfied for U p (Akt/r 1 ^ k+1 \ . . . jAot/r 1 ^^ 1 ^), where the integer r satisfies 

( eN„ k Q„. k At ' " 

r > 

and the error e satisfies 



,(2p+fc + l) 1 /(fe+D 

then \\U(t) - U p (A k t/r 1 / ( - k+1 \ . . . , Aot/r^^+^YW < e. 

Proof. Since the operators Aj are anti-Hermitian, e Ajt and e Zk * + are unitary. Errors in unitary operations are 
subadditive, so the total error is at most the number of time steps used in the simulation, r, times the error in 
simulating the commutator exponential during each time step. With r time steps, the duration of each step is 
i/ r 1 /(' c + 1 ) By Lemma 8, the error in approximating e Zkt is at most e if 



T ~ e (fc+i)/(2p) W 



t < ( ,„ i i ^M/r^n ) ^ + ln ( 2 ) 2P (N P ,kQ P ,kM) k+1 , (53) 



eN p . k Q p , k At yp+k+i 



< -, (54) 



i.e., if 



(2p + fc + l) 1 /(fc+i)r 1 /(fc+i)y - r ' 

< Mp)/(k + l) (cr\ 

e V(2p + fc + l) 1 /(fe+i) ' 



1 / eN pM Q pM At \ 2p+fc+] 
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which is equivalent to the stated condition (52) on r. 

The argument that the restriction on e in (53) implies assumptions 5 and 4 of Lemma 8 is straightforward. We first 
demonstrate that N Pik Q Ptk > 1 for U p (A k t/r 1 ^ k+1 \ . . . , Ait/r 1 ' and then use this fact to show the validity of 
assumption 4. We show this by proving loose lower bounds for N Ptk and Q Ptk arid showing that the product of these 
lower bounds yields a result greater than 1. First, we use (35) to sec that 

N p , k > 2 • 5*'- 1 7V p!fc _ 1 > (2 • 5f- 1 ) fe , (56) 

where the last inequality follows from solving the prior recursion relation using N Pt i > 2 • 5 P ~ 1 as the initial condition. 

Since Q Pjk > max.,- { | A j | } , where Xj is the duration of the j th exponential in U p (A k t/r 1 ^ k+1 \ . . . , Ait/r 1 / < - k+1 ' > ) di- 
vided by t/r l /( k+l \ a lower bound on |Aj| is also a lower bound for Q Pjfc . The form of U p (A k t j r 1 1 . . . , Ait / V 1 /( fc+1 )) 
requires that we recursively use the results of Theorem 2 and Theorem 5. Note that this recursion must be used k 
times and each such formula is constructed by p applications of the recursive relations in Theorem 2 or Theorem 5. 
Theorem 2 and Theorem 5 imply that each A.,- must be a product of at most qk terms of the form /i^, vg, fie, and 
"fi for I G {1, . . . ,p}. Because fj,g, isg, (3i, and 7f are monotonically decreasing functions of £, taking their limit as 
I -)• oo shows that max{nt,u t ,p t ,'yt} > l/3 1/(k+1) > Furthermore, £ k = 2- 1 ^ k+1 ^ > 1/-/3. Because Xj is a 

product of qk such terms, it follows that 

Q P ,k > min| A, | > S^ 2 . (57) 



Combining (56) and (57) gives 



(2 . §P~ 1 ) k 2 

N P , k Q P , k > 3kp/2 > > -f= > 1. (58) 



Since N Ptk Q p . k > 1, assumption 4 of Lemma 8 is implied by 



r i/(fc+i) - N Pik Q Ptk 
Using our lower bound on r in (52), we find that (59) is implied by 



Ate 1 /^ ln(2) , . 



( eN Pik Q Pik At \ 1+ 2 P N p , k Q p! k 
V(2p+fe+l) 1 /(fc + i) J 

which is equivalent to (53) after solving for e and simplifying the resulting expression. □ 

Rather than fixing p and choosing r such that the error is at most e, a more cost-effective strategy is to choose p 
to minimize the total number of exponentials. The required number of exponentials to achieve error at most e is 

, tin ( fc +!) 2 
/ eN P}k Q p , k At Y+ L + 2 P 

N„ p = N r . ir >N j-^'^l rt . (61) 

This expression is minimized for some p due to two competing tendencies: e~( fc +i)/(2p) shrinks as p increases, whereas 
N Pik grows exponentially in p for k constant, as discussed at the end of Section III. We can approximate the optimal 
value of p by equating two terms: 

/ M \ C=+i) 2 /(2p) 

Kt 2 - y^) ■ 

Here we neglect the term Np k k +1 ^ ^ 2p ^ since it is approximately constant for N P:k € e e ( p \ We could also include the 
detailed behavior of Q p . k to more accurately estimate the optimal value of p, but doing so complicates the discussion 
and does not qualitatively change the result. 

Taking logarithms of both sides of (62), using N Ptk € e e ( p ° pt \ and considering constant k, we find 

Popt e e L /log H . (63) 
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Since p opt is sublogarithmic in Af/e, we find that N p ^Q P ,k grows subpolynomially. Thus, by choosing p = p pt, the 
number of exponentials used to approximate e Zkt + scales as 

/AA o(1) 

^Vcxp G (At) k+1 (™) . (64) 

Similarly to the simulation of sums of Hamiltonians using Suzuki formulas [5, 6], this shows that the cost of simulat- 
ing an operator exponential scales only slightly superlinearly with the total evolution time t k+1 and subpolynomially 
with 1/e. (Note that the analogous scaling for the Suzuki formulas [ ■] follows upon substituting k = 0.) In the next 
section we show that this behavior is nearly optimal. 



VI. OPTIMALITY PROOF VIA QUANTUM SEARCH 



We have shown that high-order approximations to exponentials of commutators can be constructed with a number 
of exponentials that scales almost linearly with the evolution time. Here we show that such schemes are nearly 
optimal by demonstrating that a simulation using a sublinear number of exponentials would violate the lower bound 
for quantum search [ ] . The argument is similar to the analysis of the continuous-time algorithm for quantum search 
[20] , but uses a Hamiltonian expressed as a commutator instead of as a sum. 

Theorem 10. Any product formula approximation for e vA,B ^ T for generic A,B consists of Q(T) exponentials. 
Proof. Consider the problem of searching for an unknown w € {1, . . . , n} using a black box acting as 

, ,\ , J \x,b) if x ^ w 

\x,b)^i: ' (65) 
I \X, 0/ II X = w. 

A quantum computer must make fi(y / n) queries to the black box to solve this problem with bounded error [19]. 

To solve the search problem, we can begin in the state |+) := X)"=i \ x ) an d try to evolve with the Hamiltonian 
i[|u>)(w| , |+)(+|] for some time T such that 

e [\ W ){ W \,\+)(+\]T | +) ^ ]w) (66) 

for any w £ {1, . . . , n}. 

The Hamiltonian i[|w)(w| , |+)(+|] generates a rotation in the plane spanned by \w) and |+). We therefore choose 
our basis to contain \w) and |_L) := Y^x^w \ x )- Let 

Y:=Vn[\w)(w\,\+)(+\] 

= \w)(+\-\+)(w\. (67) 

The action of Y on the basis states is 



Y\ w ) = -d—\±) (68) 



Y\±) = J— \w). (69) 
V n 

This is simply iy/ (n — 1) /n times the Pauli operator a y on the relevant subspace, with eigenvectors 

|e ± ) = (70) 

and corresponding eigenvalues 



e± =±iJ——. (71) 
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We have 



Defining <f> so that tan</> = y/n — 1, the time-evolved state is 



s * |+> = M e *(^F*-*) 



72 



(74) 



The evolution then yields the desired state, \w), when the phases of both eigenvectors are the same. This first occurs 
when 

r = 7 ^= = fvs+0(i). (75) 

yn - 1 2 

The exponential of |u>)(u>| for any time can be simulated with 2 queries to the black box, and the exponential 
of |+)(+| for any time can be simulated without querying the black box. Thus a product formula approximation of 
e [|™X™>U+X+l] T using A^cxp exponentials gives an algorithm for the search problem using 0(N exp ) queries. By the Q(y/n) 
lower bound on search, we find that the number of exponentials in the product formula must be -/V cxp € il(y / n) = Cl(T). 

a 

Thus there is no product formula approximation to e^ 4 ' 5 !* + , for generic A and B, that consists of o((At) k+1 ) 
exponentials. Equation (64) shows that our high-order product formula approximations use (At) k+1+ °^ exponentials, 
so they are nearly optimal. 



VII. APPLICATIONS 



A. Controlling Quantum Systems 

Quantum control is the art of designing sequences of quantum operations (often referred to as pulses) that ap- 
proximate a desired operation. This is especially useful when the target operation cannot be directly implemented. 
Techniques that are used to design pulse sequences include [14] Lie-Trotter-Suzuki sequences, Dyson series, Fer and 
Wilcox expansions, Solovay-Kitaev sequences, and the Magnus series. The latter three techniques explicitly involve 
implementing exponentials of commutators, which can be turned into pulse sequences using product formulas. Cur- 
rent approaches to implement such commutator exponentials are typically inaccurate and often must be numerically 
optimized via a method such as GRAPE [21, 22], which uses gradient ascent to numerically optimize the efficiency 
and robustness of a pulse sequence that implements a desired unitary operation. The quality of the numerically 
optimized sequence depends on the initial input sequence, so our results may be useful for designing pulse sequences 
(and perhaps understanding why a numerically optimized control sequence takes a particular form). 

A simple application of our results to quantum control addresses the problem of performing an arbitrary single-qubit 
rotation given access to two non-orthogonal rotation axes. Consider the single-qubit Hamiltonian 

H(t)=B a z +B c (t)(<T z + a x ) (76) 

where B c (t) is a controllable magnetic field that can switch between the values and B c < Bq. The goal is to perform 
the operation 

U(t) = e- l " 0CT »* 2 (77) 

for some constant cjq, which is sufficient to establish universal control of the qubit if used in concert with the z rotation 
provided by the Hamiltonian when B c = 0. We can achieve the desired y rotation as 

U(t) = g- w ok z ,CT ;E ]t 2 /2 _ e -uj [Boa !: ,B i7 z +B c (a z +a- a: )]t 2 /(2B B C ) 
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For simplicity, define 

A := B Q a z , 

B:=B Q a z + ^(a z + a x ), (79) 

which can both be realized by appropriate choices of B c . Then [A, B] = iuj^Uy. Corollary 3 shows that we can 
approximate U(t) as 

U(t) = e -^t/V2 e ~iBt/V2 e iAt/y/2 e iBt/V2 e iAt/V2 e iBt/V2 e -iAt/V2 e ~iBt/^ + (J(t A ). (80) 

Theorem 2 can be used to construct higher-order variants of this expression. 

Note that we have assumed that we can implement evolutions of the form exp(iAt) and exp(iBt) for t > 0, 
corresponding to evolution for negative time. For this single-qubit system, this can easily be done by evolving under 
A or B for some positive time. However, this may not be as straightforward in general. 

Recent work by Borneman, Granade, and Cory [23] further highlights the role of commutator exponentials in 
quantum control by using them to suppress undesirable interactions in coupled spin systems and introduce new 
couplings that are not present in the original Hamiltonian. Their work constructs a nested commutator of the form 
e t[A,[A,B]\t usm g on iy the ability to implement exponentials of A and B individually. They achieve this by presenting 
a product formula for nested commutators that is outside the scope of our formalism (since it exploits the fact that 
two of the terms in the nested commutator are equal) , namely 

W$$(-iAt, -iBt) := e -iAt e -iBt e iAt e iBt e iAt e -iBt e -iAt e iBt_ (81) 

They show that this product formula obeys 

Wf% c 2 {-iAt, -iBt) - e ^.[^]]t 3 || e ( t 4), (82 ) 

Here we use the subscript p — 1/2 since, with k = 2, the error is 0(t 4 ) = 0(t 2p+k+1 ) as in our previous product 
formulas. The validity of (82) follows directly from 

ln(e- iAt e- iBt e iM e iBt ) = -[A, B]t 2 + ^[A + B, [A, B}} + ^([A, [B, [B, A}}}/2 + [A + B, [A + B, [A, B}}]) + 0(t 5 ) 

(83) 

and the Baker-Campbell-Hausdorff formula because all terms of 0(t 3 ) cancel in (81) except for it 3 [A, [A, B]]. 

Note that W^B^(—iAt, —iBt) is actually correct to 0(t 5 ) if [A, [B, [B, A}}} — 0, which holds, for example, when A 

and B are Pauli operators. We refer to the BGC formula as W®2 C {~ iAt, —iBt) in such cases to reflect the fact that 
the error is 0(t 5 ). Also note that by a similar calculation to (83), 

Wfj$(-iAt, -iBt)' 1 = e -iBt e iAt e iBt e -iAt e -iBt e -iAt e iBt e iAt (g4) 

satisfies 

\\W$$(-iAt, -iBty 1 - e -i\MA,B)]t 3 \\ G (t 4 ). (85) 

This error term is also 0{t 5 ) when [A, [B, [A, B]]} = 0. 

Even though the BGC formula goes beyond our formalism, we can still use our results to refine it into higher- 
order versions. It is not immediately clear how to do this using Theorem 5 since the proof of that theorem relies on 
determining the parity of the number of As in a particular error term. For the error term of order t 5 in the BGC 
formula, we cannot determine the parity of the number of As because both A and B are associated with t rather than 
with t and t 2 , respectively. 

Nevertheless, it turns out that the construction described in Theorem 5 refines W^S g(— iAt, —iBt) into a formula 

denoted Wf$°{-iAt, -iBt) that has error 0(t 5 ) (rather than 0(t 6 ) for the reasons described above). We can then 
continue to refine the formula to arbitrarily high order by applying the construction of Theorem 5 recursively. In 
particular, we have the following. 

Corollary 11. Let Wp^ c (At, Bt) be a product formula approximation to exp([A, [A, B]]t 3 ) that has error 0(t 2p+3 ). 
Define 

K+i%^At, Bt) := W^ c (Av p t, B Vp t) 2 W^ c (A^ p t, Bti P t)- x Wf$ c {Av p t, Bv p t) 2 (86) 
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where 

3 

Mp := (4 Sp ) 1/3 , v p := (1/4 + sp) 1 / 3 , Sp := 4>+ ' ■ (87) 

4(4 — 4 2 p+ 3 ) 

T/ien ||W; B + G ^ 2 , 2 (Ai,St) - e ^^ A ^]] t3 || S 0(t 2p+i ). 

Proof. We have % B + G 1 / 2 . 2 (^, Bt) = exp([A, [A, B]]t 3 ) + C(A, B)t 2 P+ 3 for some operator C(A, B). Similarly, from (85), 
W p+i%.2( At > Bt )~ l = cx pHA [A B ]]t 3 ) + C{A, B)t 2p+3 for some operator C(A, B). Since 

= 1 + (C(A, B) + C(A, B))t 2p+3 + 0(t 2p+i ), (88) 

C(A, B) — —C(A, B). Then the result follows from the same calculation as in (27). □ 

An example of this generalized BGC formula with A = ia x and B = i<r z is presented in Figure 3. Here the 
lowest-order formula is W /B 2 3C (— ia x t, —ia z t) since [A, [B, [B, A])] = 0. We see that the upper bound of Corollary 11 is 
essentially optimal: the approximation error scales as 0(t 2p+3 ). Note that our higher-order BGC formulas use 0(5 2p ) 
exponentials, an improvement over the 0(6 2p ) exponentials appearing in the results of Lemma 6. 

As another simple application, we could also use the nested commutator formulas provided by Lemma 6 or Corol- 
lary 11 to implement e[ j4 '^' S H t in the simple example discussed at the beginning of this section, where A and B are 
given by (79). This produces an effective x rotation, effectively canceling the a z terms present in B. 



B. Product Formula Approximations for Anticommutators 



Another application of our results is to implement Hamiltonian evolution according to anticommutators. The 
anticommutator of two operators A and B is 



{A,B} := AB + BA. 



(89) 
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We can implement e^'^* using a commutator simulation by enlarging the Hilbert space. This method can also 
be used to implement powers of a Hamiltonian, or more generally, products of powers of commuting Hamiltonians, 
which may be applied to simulate many-body couplings. In contrast to the techniques described in the previous 
section, this method can be used to introduce interactions that do not arise naturally as commutators of terms 
already appearing in the Hamiltonian. Our product formula approach for simulating anticommutators can therefore 
be seen as a complementary method to the generalization of the BGC approach described in Section VII A. 

We extend the Hilbert space from T-L to %(g>C 2 and create analogs of our operators A and B on this enlarged space, 
namely 



A 1 := A ® a y 

B':=B®a x , (90) 



where a x and a y are Pauli operators. Then 



[A 1 , B'] = AB ® a y a x - BA ® a x a y 

= -i{A 1 B}®a z . (91) 

Let G H be any quantum state. Then 

[A', B'\ |V) ® |0) = -i{A, B} |V> ® |0) , (92) 
so the action of the anticomutator can be simulated using an ancilla qubit in the eigenstate |0) of a z . In particular, 

e ^{A,B}t- |^ g | Q) = e [A',B']t* g | 0) 

= V^CA't, S't) |V) ® |0> + <3(t 2p+2 ). (93) 

This shows how to simulate anticommutators using a simulation of commutators. It is straightforward to generalize 
this construction to simulate nested anticommutators in terms of nested commutators. 

While we have already discussed applications of implementing exponentials of commutators, exponentials of anti- 
commutators may seem less natural. However, one application of anticommutators (between commuting operators) 

is to simulate many-body couplings. For example, consider the one-body operators Aj = a% for j = 0, . . . , k, where 
the superscript indicates which spin is acted on. Then 

{A k ,... {A 2l {A^Ao}} ...} = 2 k af k+ \ (94) 

so the Hamiltonian af k+1 can be simulated using a product of exponentials that each involve only two-body terms, 
which are typically more natural than fc-body terms with k > 2. This approach provides an alternative to direct 
simulation methods for many-body terms [24-27] . Since no explicit changes of basis are performed, this simulation 
may more accurately reproduce features of the ideal evolution (e.g., its behavior in the presence of errors). 
As an example of a complete simulation, consider the toric code Hamiltonian [28] 

H = -j[Y d Av + Y, B P J ' ( 95 ) 



where A v := Yiiev a ^ anc ^ := T\jep a * are four-body coupling terms acting on qubits on the edges of a square 
lattice, where v denotes the set of qubits surrounding a vertex of the lattice and p denotes the set of qubits on the 
edges of a plaquette. The vertex and plaquette operators commute, so the time-evolution operator is 



-iHt 



II' " II 



ij B n t 



(96) 



Each of these exponentials involves a four-body Hamiltonian and therefore can be implemented using only two-body 
interactions via the approach outlined in (94) with k = 3. 

The following theorem generalizes the above examples to simulate a Hamiltonian that is a product of powers 
of commuting matrices and provides a loose upper bound for the scaling of the number of exponentials (which is 
proportional to the number of quantum operations) needed to simulate the Hamiltonian. 
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Theorem 12. Let {Af. £ = l,...,m} be commuting Hermitian matrices, let {cti: I — l,...,m} be integers, let 
fc := Q^aLi a <0 — 1- an d let p be a positive integer. Then exp(— i2 k A" 1 . . . A^™t kJrV ) can be simulated with error at 
most e > using 



fe+l + (fe+l) 2 /(2p) 



N exp EO\ l | ( «)7) 



exponentials, where A > 2maxj ||j4j| 



Proof. The proof is a simple generalization of the previous discussion. First introduce operators Aj for j = 0, . . . , k 
such that each Aj corresponds to one of the operators At, where Aj — Ag for on values of j. Then define A'j as the 
isometric extension 

A , := U®* X ifi = 

J ® 0~y if J > 0. 

We now show by induction that the nested commutator of all the A'j is an isometric extension of .4o-4i . . . Ak ■ 
Specifically, we claim that 



[A' k ,[...,[A'i, ■<]■■■}] 



-*2 fe U q =i A ® cr x k odd 
2 fc A q <S> ct z fc even 

-i2 fe Ag°A" 1 • • • A£ fe ® da; fcodd 



(100) 



\2 k A%°A^ ■■■A k ' k ®cr z k even. ^ 

We have already demonstrated the base case in (92) by showing that [^4'i,^4q] = — 2iAiAo <8) a z . Thus, assume that 
the claim holds for some given value of fc. Then we have 

[A k+1 , [A k ,[..., [A^Aol . . .jjj - | 2fe+1 ^ g ^ ^ fc ^ 

_ f 2 fc+1 A ® ^ fc + 1 even 
~ \ -i2 k+l \{ ] q t\ A q ®a x fc + 1 odd, 

and the claim follows by induction. 

Finally, we simulate the nested commutator exponential by U p (A k t, . . . ,A' Q t) (with an eigenstate of a z or a x , as 
appropriate, in the ancilla register). The theorem follows by using (61) for N p . k G 0(6 pfe ) and Q Ptk < 1 (as justified 
in Section V). □ 

In the example of the toric code, Theorem 12 implies that e~ lHt can be simulated with error at most e using 

\ {e/nf/P ) {WL) 

operations, where n is the number of vertices (or equivalently, the number of plaquettes) in the lattice. We replace e 
with e/n because simulation errors are subadditive and there are 0(n) terms to be simulated in (96). Equation (101) 
shows that the toric code dynamics can be implemented using a number of two-body interactions that scales near- 
linearly with t while using only one ancilla qubit (although it may be more convenient to use 0(n) ancilla qubits so 
that nearest-neighbor interactions suffice). 

An alternative to Theorem 12 is to estimate the eigenvalues of each Ai using phase estimation and introduce a 
conditional phase that depends on the eigenvalue (e.g., as described in [29]). The method based on commutators has 
the advantage that it needs only one ancilla qubit, rather than a logarithmic-size register used to store the estimated 
phase. Also, phase estimation algorithms have the drawback of requiring 0(l/e) operations in general. In contrast, 
our approach uses e -0 ^ 1 ) operations, so Theorem 12 may be useful for high-precision simulations. 
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VIII. CONCLUSION 

We have presented recursive constructions that approximate exponentials of commutators as products of exponen- 
tials of the elementary terms. We provided explicit upper bounds on the approximation error, found upper bounds on 
the number of elementary exponentials needed to approximate an exponential of commutators to within a fixed error 
tolerance, and established near-optimality of our results by relating commutator simulation to unstructured search. 

Our formulas have natural applications in quantum control, where they can be used to introduce or suppress 
interactions. Our formulas can be much more accurate than those used in current applications, and consequently may 
have applications for designing highly accurate control sequences or dynamical decoupling sequences. Such sequences 
may be valuable as initial guesses in numerical pulse finding methods. 

Anticommutators can also be simulated using commutator simulation together with a dilation of the Hilbert space 
to switch the signs of the terms in the anticommutator. This technique can be used to implement many-body 
Hamiltonians or simulate powers of Hamiltonians. 

This work raises several natural questions. Although our product formula approximations have nearly optimal 
scaling with t, alternative methods could come closer to or even saturate the lower bound. However, this problem 
may be difficult as it has long been open in the case of exponentials of sums. A more tractable goal might be to 
seek improved performance as a function of k, which would be useful since the number of exponentials needed for a 
high-order approximation for large k can be prohibitively expensive. 

Applications of our results to quantum control are reminiscent of the Solovay-Kitaev theorem [30], which is also used 
to generate pulse sequences. That approach is based on the group commutator, which we also use in Lemma 1. Our 
expression in (80) has the advantage of more precisely implementing the desired rotation than the group commutators 
used in Solovay-Kitaev decompositions [31-33]. This suggests that our product formulas may also find application in 
Solovay-Kitaev algorithms. 

Finally, recent work in numerical analysis [34, 35] and quantum simulation [36] has shown that multi-product 
formulas, which are linear combinations of product formulas, can provide more efficient approximations to operator 
exponentials. It might be interesting to investigate whether similar techniques could lead to more efficient (but 
possibly less numerically stable) approximations to exponentials of commutators. 
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